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Abstract 

Temporal evolutions toward thermal equilibria are numerically investi- 
gated in a Hamiltonian system with many degrees of freedom which has 
second order phase transition. Relaxation processes are studied through lo- 
cal order parameter, and slow relaxations of power type are observed at the 
critical energy of phase transition for some initial conditions. Numerical 
results are compared with results of a phenomenological theory of statis- 
tical mechanics. At the critical energy, the maximum Lyapunov exponent 
takes the largest value. Temporal evolutions and probability distributions 
of local Lyapunov exponents show that the system is highly chaotic rather 
than weakly chaotic at the critical energy. Consequently theories for per- 
turbed systems may not be applied to the system at the critical energy in 
order to explain the slow relaxation of power type. 

1 Introduction 



In recent years some researchers have investigated Hamiltonian systems with 
many degrees of freedom |l| (and references therein). They are interested in 
statistical mechanics since Hamiltonian dynamics is the basis of it. They mainly 
have considered relaxation processes to equilibrium states because we can observe 
temporal evolutions of systems directly from Hamiltonian dynamics. 

In both the one-dimensional classical $ 4 lattice model and FPU (3 model, 
they observed slow and fast relaxation processes in low and high energy regions, 
respectively The slow and fast relaxation processes occur in weakly and 

highly chaotic systems, respectively M. Here the word of weakly chaotic system 
means that theories for perturbed systems may be applied for the system. For 
instance, if we add small non-integrable perturbation to a completely integrable 
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system, whose phase space consist of tori, almost all tori survive under the small 
perturbation. That is a result of Kolmogorov- Arnold- Moser (KAM) theorem 0]. 
Intermittency between chaotic and regular-like motions is also a characteristic 
phenomenon which occurs in perturbed systems. Arnold diffusion ,|| Nekhoro- 
shev timeJTJ and other theories exist for perturbed systems. Using theories and 
properties of weakly chaotic systems, above-mentioned slow relaxation, which 
occurs in weakly chaotic systems, was discussed^. 

Properties of weakly chaotic systems have been well investigated both theoret- 
ically and numerically. Contrary to weakly chaotic systems, few studies exist for 
highly chaotic systems, and to study properties of highly chaotic systems whose 
dynamics can not reduced to Markovian process is a left important problem of 
classical Hamiltonian systems. 

We expect it is useful to consider relaxation processes in a Hamiltonian system 
having second order phase transition in order to study properties of highly chaotic 
systems. The reason is as follows. The system seems highly chaotic near the 
critical point because fluctuation anomalously increases. Nevertheless, according 
to a phenomenological theory of statistical mechanics, van Hove theory, slow 
relaxation of power type of order parameter appears at the critical point. That 
means a certain structure exists in phase space. 

In this way we understand that critical phenomenon is interesting as a problem 
in dynamics, as well as in statistical mechanics. 

In this paper, we investigate the following things. First, as a problem of sta- 
tistical mechanics, we verify appearance of the slow relaxation of order parameter 
at the critical point even when we use a dynamical method, and we compare re- 
sults of dynamics with ones of phenomenological theory. Second, as a problem of 
dynamics, we confirm that the system is highly chaotic at the critical point, and 
we suggest that critical phenomenon is a useful phenomenon to study properties 
of highly chaotic systems. 

We introduce a Hamiltonian system in which second order phase transition oc- 
curs, and numerically integrate the equations of motion. The model Hamiltonian 
is as follows || 

N i 

%p) = E-i + p(!). (i) 
»=i A 

i N 

U(q) = ^rJ2( 1 - C0S (li-q j ))- (2) 

Equations of motion are derived from the Hamiltonian as follows 

dqi , * 

dp- 1 N 

dT = -ivE^-g,), (4) 
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and temporal evolutions of the system are yielded from the equations of motion. 
Each particle moves on the unit circle and interacts with all the others. The g^'s 
are phase of particles and the p^s are canonical conjugated momenta. Accord- 
ing to statistical mechanics, the critical energy E c is E c /N = 0.75[Q, where N 
represents degrees of freedom. 

We define an order parameter parameter M of this system as 

M= \ [ L dtMCt) (5) 
L Jo 



and 

M{t) = \\M{t)\\, M(t) 



' 1 N 1 N \ 

— 5>osg;(t), — J2 sin Qi(t)J > ( 6 ) 



where t and | \M(t) 1 1 represent time and absolute value of M(t), respectively. From 
the definition, < M, M(t) < 1. When all particles distribute uniformly on the 
unit circle at each time, then M = 0. On the other hand, when particles form the 
cluster, which is the only one in this system |L0| , then M ~ 1. For this system, 
Antoni and Ruffo showed that energy dependence of order parameter obtained 
from numerical experiments agrees with a result of statistical mechanics || . 

We observe relaxation processes, that are temporal evolutions of M(t), toward 
equilibria using a local order parameter M T (n), which is defined as 

1 r nT 

M T {n) = - / dtMCt). (7) 

T J(n-V)T 

The system Eq.(|l|) is integrable in the limit E/N —>■ and E/N — > oo. For 
the former and the latter limit the system is approximated as harmonic oscillators 
and free particles, respectively. To explain the case of the latter limit, we show 
the following expression of the potential term U(q): 

N 1 N 

U{q) = — - — (cos & cos g,- + sing i sing i ) (8) 

= ^(l-M(t) 2 ), 0<M(t)<l (9) 

(cf. Eq.@). Hence U(q) has the minimum and the maximum value 

< U(q)/N < 1/2. (10) 

We can therefore neglect the potential term when the value of E/N is large 
enough, then the system consists of free particles. Note that Ai —>■ when 
E/N -> and E/N -> oo (see inset of Fig.|). 

To compare our numerical results with results of a phenomenological theory of 
statistical mechanics, we introduce a phenomenological theory: van Hove theory. 
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We assume temporal evolution of order parameter M(t) is determined by gradient 
of Landau's free energy F(M), that is, 

dM(t) dF(M) 

-dr"—djir> 

F(M) = a'(T -T C )M 2 + b(T)M A , (a' > 0). (12) 
Note that relaxations to equilibria occur at once if states are non-equilibria 



(dF(M)/dM ^ 0). From the solution of Eq.flTTj), we obtain the following tem- 
poral evolution of M{t) 

/ eX P( _t / T ( T )) ( When T < T c) f.n\ 

M (*)~{ rV2 (when T = T C ) ' (13) 

where t(T) is called relaxation time which depends on temperature T and di- 
verges at critical temperature T c , and Eq.([H|) is approximately correct for t 3> 
r(T) when T < T c . According to this theory, relaxation of power type appears if 
and only if systems are just on critical points. 

To integrate the equations of motion which are derived from the Hamiltonian 
of Eq.(j|D, we use forth order symplectic integrators [|ll|] with fixed time slice 
At = 0.01 (At = 0.001 for a part of result in Fig.^) which keep symplectic 
properties of Hamiltonian systems exactly and total energy accurately. Relative 
errors of total energy AE/E are less than O(10 -7 ) for considered energy region. 
The value of r in Eq.(^) is fixed at 10, namely 1000 steps except for the case 
given special comments. As initial conditions we choose random variables which 
follow Gaussian distribution. To see relaxation of order parameter, we take small 
values for q^s and set M ~ 1. Scales of pj's are defined from energy. In the 
model Eq. (H) total momentum is conserved, and we set total momentum is equal 
to zero. 

This paper is organized as follows. In Section 2 results of numerical exper- 
iments are reported. At first, slow relaxation of power type is shown, then we 
investigate dependences on degrees of freedom and initial condition for appear- 
ance of the slow relaxation. After that we compare our numerical results with 
results of van Hove theory. Next we confirm that the system is highly chaotic at 
the critical energy with the aid of temporal evolutions and probability distribu- 
tions of local order parameter. Section 3 is devoted to summary and discussions. 

2 Results of simulations 

In this section numerical results are reported for the Hamiltonian system Eq. ([!]). 



In Sec.|2.1| we find that slow relaxation of power type appears at the critical energy 



for a certain degrees of freedom and initial condition. We observe dependences on 
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degrees of freedom and initial condition for appearance of the slow relaxation in 



Sec.|2~^ and Sec|Q|, respectively. Section 2.4 is devoted to compare our dynamical 
results and results of a phenomenological theory of statistical mechanics. After 
that we investigate dynamical properties of the system, in particular, at the 
critical energy. Energy dependence of maximum Lyapunov exponent is shown in 



Sec. |2.5| . Then, in Sec. [2.6| , we confirm that the system is highly chaotic at the 
critical energy. That is, the slow relaxation of power type may not be explained 
by theories for perturbed systems. We study the structure of phase space in 
Sec.|2~?l, and suggest phase space is uniform at the critical energy for the property 
of instability. 



2.1 Slow relaxation at the critical energy 

We show a energy dependence of order parameter M in Fig.|l|. The solid line 
in the figure represents the curve obtained from theory of statistical mechanics. 
Using the saddle point method, it is described as simultaneous implicit functions 

= 

M=^(J3M), (14) 

E/N=^+ l -{l-M% (15) 

where Jo and I\ are Bessel functions of 0-th and 1-st order, respectively. Then we 
find that results of numerical experiments fit the theoretical curve, and the results 
are reproduction of results of Antoni and Ruffo||. Although we find differences 
between the results and theory, we can understand causes of the differences as 
follows. For high energy part (E/N > 1), M(t) fluctuates around zero, which is 
the minimum value of M(t), since the system is approximated as free particles, 
and <7i(i)'s takes random variables. Then M(t) is estimated around 0{l/yN) 
from central limit theorem. That is confirmed from N dependence of the values 
of M in this part. For middle energy part (E ~ E c ), in addition to the cause 
of differences in high energy part, there is another cause. That is lack of time L 
in which we take time average of M(t) because, according to van Hove theory, 
relaxation time increases as energy goes to critical value. 

Here let us investigate relaxation of M T {n) defined as Eq. ([?[). We report results 
when E/N = 0.5, 0.75(= E c /N) and 1 in Fig.^. The insets of Fig.|2] represent the 
following quantity as a function of time t 

M(t;t ) = — !— fdtM{t). (16) 
t — to Jt 

This quantity represents time average of M(t) from to to t, and reduces fluctuation 
£(£) when temporal evolution of M(t) is as follows 

M(t) = (t~t )~ x + at)- (17) 
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For Fig.|2](b) we set to = 1000 since slow relaxation of M T (n) starts around the 
to- For Figs.^(a) and (c) we set to = 0. From Fig.^, we confirm the relaxation 
is power type when E = E c for the initial condition. This result agrees with a 
result of van Hove theory. The slow relaxation finishes in finite time since M(t) 
takes degree of 0(1/ VN) for t — > oo. 

2.2 N dependence of slow relaxation 

We observed slow relaxation of power type when N = 80. Now we show that 
the slow relaxation appears even if we change degrees of freedom. When N = 40 
the slow relaxation of power type is shown in Fig|| Inset of the figure represents 
M(t;io) (cf. Eq.(|n5D against t again, where we set to = 5000. However we have 
not observed the slow relaxation when N > 80, and we will discuss the reason 
and behaviors of the system in thermodynamic limit in Sec.|3| 

2.3 Initial condition dependence 

In the previous section, we investigated dependence on degrees of freedom N. 
Since we treat dynamical system we must study dependence on initial conditions 
for relaxation of M T (n), too. A problem is whether slow relaxation always appears 
or not when E = E c . Figure |] shows results when N = 80, 200 and 1000. 
Exponential relaxations become clearer and clearer as N increases hence slow 
relaxation of power type does not always appear. A rate of initial conditions 
which yield slow relaxation may be low since only one initial condition yields 
slow relaxation in ten initial conditions when E = E c and iV = 80. 

2.4 Comparing with a phenomenological theory 

Up to now we observed appearance of slow relaxation of power type at the critical 
energy for some initial conditions although the slow relaxation does not always 
appears for any initial conditions even if E = E c . The appearance of the slow 
relaxation is a agreement with a result of a phenomenological theory of statis- 
tical mechanics: van Hove theory. However a disagreement also exists between 
our results and results of the phenomenological theory, and here we discuss the 
disagreement. 

We found there are two time regions for slow relaxation process in Figs.|2|(b) 
and One is the region in which the system stays in non-equilibrium states (flat 
region), the other is the region in which slow relaxation of power type occurs 
(relaxation region). The existence of the flat region is different from theory of 
van Hove, which says systems in non-equilibrium states goes toward equilibrium 
states immediately (cf. Eq. flTTD ). Our dynamical results are therefore an example 
of a disagreement with a result of van Hove theory. 
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The flat region means the existence of "induction period" and that reminds 
induction phenomenon, which means that a nearly periodic motion starts to 
behave as a non-periodic motion after a certain amount of time, which is called 
induction period |T2[||13||P^ |. However the phenomenon appearing in Figsf2|(b) 



and |3] is different from induction phenomenon since the motion in phase space is 
not nearly periodic both before and after the onset of slow relaxation. To confirm 
that we show average of power spectra of momenta pj's (j = 1, 2, • • • , No), S p (f), 
when iV = 80 in Fig.|. The quantity S p (f) is defined as 



I N 



W) = TrEW), (18) 

^0 3=1 



Sj(f) = "power spectrum of time series of Pj(t)" . (19) 

Here we selected iVo momenta from N momenta to calculate S p (f), and we set 
iVo = 40. In low frequency region the power spectrum is growing up hence the 
motion is not nearly periodic. 

2.5 Lyapunov exponent 

We have investigated slow relaxation of power type using dynamical method. 
Now we study dynamical properties of the slow relaxation and structure of phase 
space in which the slow relaxation occurs. 

As the first step, we investigate an energy dependence of the maximum Lya- 
punov exponent Ai (hereafter we call this Lyapunov exponent simply) which 
measures instability of orbits. Existence of positive Lyapunov exponent means 
a sample orbit has instability. Results of numerical experiments are reported in 
Fig.|| and we find Lyapunov exponent takes the largest value at E c /N{= 0.75) 
for all sample orbits when N = 40, 80 and 200. 

We can understand why value of Lyapunov exponent takes the maximum at 
E c as follows. The system is integrable for low and high energy limit since it is 
harmonic oscillators and free particles, respectively. Thus Ai —>■ when E/N — ► 
and E/N — ■> oo. Integrabilities of those two limits break as energy goes away 
from zero or infinity, and the two integrabilities balance at the critical point. 
Consequently Lyapunov exponent is the maximum at the critical point. 

Butera et al. also investigated Lyapunov exponent for a Hamiltonian sys- 



tem which is a two-dimensional system of coupled rotators [15]. They found a 
"knee-like" shape in the graph of Lyapunov exponent against temperature and 
the "knee" is at the critical temperature of Kosterlitz-Thouless (KT) phase tran- 



sition |16| although the "knee" does not appear at the largest value of Lyapunov 
exponent. Figure ^| shows a "knee-like" shape and the "knee" is at the critical 
energy where Ai is the largest. This fact suggests that degree of instability has 
some relations to phase transition. 
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2.6 Highly chaotic property at the critical energy 

We found that Lyapunov exponent takes the largest value at the critical energy 
hence we expect the system is highly chaotic at the critical energy. In this section 
we confirm that the system is highly chaotic at the critical energy from tempo- 
ral evolutions and probability distributions of local Lyapunov exponents. The 
definition of the local Lyapunov exponent is as follows 

\ T (n) = - / dt\i(t), (20) 

Ai(f) = ^log(||X(t)||) (21) 

where X(t) is a 2iV-dimensional tangent vector at time t which obeys linearized 
equations of motion. The local Lyapunov exponent X T (n) indicates instability of 
orbits in the time interval [(n — l)r, nr]. Orbits are instable if X T (n) > in the 
time interval. 

Figure [7| shows temporal evolutions of local Lyapunov exponents for E/N = 
0.2 and 0.75(= E c /N). The initial condition of Fig.[7|(b) yields slow relaxation of 
power type (cf. Fig.||(b)). When E/N = 0.2, typical intermittency is found. On 
the other hand, when E/N = 0.75 = E c /N, behavior of the temporal evolution 
is different from the case of E/N = 0.2 hence motion is not intermittency for 
the initial condition. Moreover local Lyapunov exponents do not take near zero, 
hence few tori, if any, exist in phase space. Consequently the system is highly 
chaotic. 

To make sure of the difference between the two time series of local Lya- 
punov exponent, we show probability distributions of local Lyapunov exponents 
in Fig|| If the system is weakly chaotic and intermittency occurs, a probability 
distribution of local Lyapunov exponents should have a peak near zero. However, 
according to Fig.[8], the distributions have not peaks near zero at E — E c . Thus 
we understand the system is different from weakly chaotic system at E = E c . 
Furthermore, at E = E c , motion is not reduced to Markovian process since re- 
laxation of power type is observed. 

2.7 Uniformity of phase space at the critical energy 

We observed slow relaxation of power type in highly chaotic system. If the system 
is weakly chaotic, we can understand the cause of the slow relaxation from the 
structure of phase space using theories for perturbed systems. However the slow 
relaxation occurs in highly chaotic system, thus we can not use the theories and 
we must newly investigate the structure of phase space to understand the cause 
of the slow relaxation. Here we suggest that phase space is uniform at the critical 
energy. 
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For the purpose, we calculate Lyapunov exponent for ten initial conditions 
when E = E c and N = 80. Results follows; 



That is, 



Ai = 0.229331, 0.229133, 0.231313, 0.229381, 0.229121, 

0.231387, 0.230692, 0.231396, 0.228699, 0.230357. (22) 

Ai = 0.2300 ±0.0014 (E/N = 0.75, iV = 80). (23) 



N standard deviation 



40 
80 
200 



0.065916 
0.056067 
0.055623 



On the other hand, Table [I] shows stan- 
dard deviations of local Lyapunov ex- 
ponents for the critical energy E c /N = 
0.75 described in Fig|| Then we find 
that the distribution of Lyapunov expo- 
nents is narrow enough comparing with 
the standard deviations of local Lya- 
punov exponents. Hence initial condi- Table 1: Standard deviations of lo- 
tions have no influences on the values of cal Lyapunov exponents when E/N = 
Ai, in other words, phase space is uni- EJN = 0.75. Those probability disci- 
form for the intensity of instability at butions are described in Fig.|. 
the critical energy. 

The values of Ai are not different between relaxations of power type and 
exponential type because the first and the second value of Eq.(p2|) correspond 
to the initial conditions of Figs.^|(b) and f|(a), respectively. Hence we can not 
distinguish between the two types only from the values of Ai, which are time 
average of instability. 



3 Summary and discussions 

Relaxation processes are numerically investigated through temporal evolutions 
of local order parameter in a Hamiltonian system which has second order phase 
transition. Slow relaxations of power type of order parameter are observed at 
the critical point for some initial conditions. It is the first time that power 
type decay is observed for order parameter, which is a important quantity to 
observe the system, from temporal evolutions of equations of motion derived 
from a Hamiltonian. 

To understand the cause of the slow relaxation dynamically, we investigated 
dynamical properties of the system. Then we found that the slow relaxations 
occur in highly chaotic systems rather than weakly chaotic systems. That is, 
mechanism of the slow relaxations can not be explained by theories for perturbed 
systems, hence we must consider new theories for highly chaotic systems. 
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The slow relaxations of power type are observed when N = 40 and 80. How- 
ever we have not observed the slow relaxation when N > 80, and some initial 
conditions yield exponential relaxations even if E = E c . 

Why the slow relaxation has not been observed when N is large? Large 
number of N makes a situation such that to detect slow relaxation is difficult 
because the larger N is, the clearer the critical point may be, and energy has 
small errors in calculations although we set E = E c initially. Moreover we do 
not know the exact value of the critical energy in the meaning of dynamics. In 
other words, we do not know whether dynamics also has the same critical energy 
E c which is obtained from statistical mechanics. However to obtain exact value 
of critical energy is difficult since appearance of the slow relaxation depends 
on initial conditions. We can not say whether the slow relaxation appears in 
the thermodynamic limit. The author expects the appearance because we have 
supporting evidences for appearance but do not have for absence. 

Appearance of the slow relaxation agrees with a result of phenomenological 
theory of statistical mechanics; van Hove theory. However, contrary to the phe- 
nomenological theory, relaxation processes do not start at once even the system 
are in non-equilibrium states. That is, we observed an agreement and a dis- 
agreement between our results from numerical experiments and results from a 
phenomenological theory. To research the origin of the agreement is interesting 
to understand contracted dynamics which yields motion of order parameter. 

From temporal evolutions and probability distributions of local Lyapunov ex- 
ponents, we understood the system is highly chaotic rather than weakly chaotic 
at the critical energy. Hence we may not apply theories for perturbed systems 
to understand the slow relaxation of power type at the critical energy. We must 
make new theories for highly chaotic systems. For the purpose we must under- 
stand structure of phase space when the system is highly chaotic and has many 
degrees of freedom. To understand the structure of phase space, a result is ob- 
tained such that the structure of phase space may be uniform for the intensity of 
instability at the critical energy. 
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Figure 1: Order parameter M vs. energy per unit particle E/N. L = 2 23 x 0.01 
(~ 10 5 , namely 10 7 steps) (cf. Eq.(|)). «»: N = 40. (+): N = 80. (□): 
N = 200. Solid line represents a result of theory of statistical mechanics. Results 
are good agreement with the theory. 
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Figure 2: Log-log plotted temporal evolutions of M T {n) (cf. Eq.(|7])). iV = 
80, t = 10. (a) E/N = 0.5. (b) E/N = E c /N = 0.75. (c) E/N = 1. Insets: The 
temporal evolution of M(t;to) which is yielded by the time series of M T (n) (cf. 
Eq.(|T6])), where t = 1000 for (b) and t = for (a) and (c). When E = E c slow 
relaxation of power type appears from t . 
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Figure 3: The same with |(b), but N = 40, t Q = 5000. 
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Figure 4: Log-log plotted temporal evolutions of M T {n) (cf. Eq.(|7|)) at E c . r = 10. 
Insets: Semi-log plotted graphs, (a) iV = 80. The initial condition is different 
from one of Fig.|(b). (b) N = 200. (c) N = 1000. We find relaxations are 
exponential type rather than power type. 
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Figure 5: Average of power spectra of momenta pj's (j = 1,2, •••,40) when 
N = 80 and the initial condition yields slow relaxation of power type. Lower 
and upper graphs are calculated from time series before (t G [0, 1024]) and after 
(t G [2048,3072]) the onset of slow relaxation, respectively. Distinguishing the 
two graphs, values of the vertical axis are multiplied 10 times for the upper graph. 
The motion in phase space is not nearly periodic since the graphs do not consist 
of solitary peaks. 
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Figure 6: The maximum Lyapunov exponent Ai vs. energy per unit particle E/N. 
L = 2 23 x 0.01 (~ 10 5 , namely 10 7 steps) (cf. Eq.(|)). Inset: The horizontal axis 
is logarithm scale, r = 0.001, L = 2 20 x 0.001 when E/N > 10 3 . (0): N = 40. 
(+): iV = 80. (□): N = 200. The quantity Ai takes the largest values at critical 
energy, namely E c /N = 0.75, and it goes to zeros in the limit of E/N — > and 
E/N oo. 
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Figure 7: Temporal evolutions of local Lyapunov exponents. N = 80. (a): 
E/N = 0.2, r = 50. (b): E/N = 0.75 = E c /N, r = 10. The initial condition of 
(b) yields relaxation of power type. Typical intermittency is found in (a), and 
(b) is not intermittency. 
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Figure 8: Probability distributions of local Lyapunov exponents. (^>): N = 40. 
(+): N = 80. (□): N = 200. The value of energy is E/N = 0.75 = E c /N. 
Dashed line is the case iV = 200, E/N = 0.2, and intermittency occurs. We set 
r = 100 when E/N = 0.2. At the critical energy the system is different from 
intermittency. 
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